Kinetic theory of nonequilibrium stochastic long-range 
systems: Phase transition and bistabihty 



Cesare Nardini^'^, Shamik Gupta\ Stefano RufFo^'^, Thierry 
Dauxois^ and Freddy Bouchet^ 

^ Laboratoire de Physique de I'Ecole Normale Superieure de Lyon, Universite de Lyon and 

CNRS, 46, Alice d'ltalic, F-69007 Lyon, France 

^ Dipartimento di Fisica e Astronomia and CSDC, Universita di Firenze and INFN, via G. 
Sansone 1, 1-50019 Sesto Fiorentino (FI), Italy 

^ Dipartimento di Energctica "Sergio Stecco" and CSDC, Universita di Firenze, CNISM 
and INFN, via S. Marta 3, 50139 Firenze, Italy 

E-mail: cesare .nardiniSgmail . com, shamikglOgmail . com, 

stef ELno . ruf f oSgmail . com, thierry . dauxoisQens-lyon . f r , f reddy . bouchetOens-lyon . f 

Abstract. We study long-range interacting systems driven by external stochastic forces 
that act collectively on all the particles constituting the system. Such a scenario is frequently 
encountered in the context of plasmas, self-gravitating systems, two-dimensional turbulence, 
and also in a broad class of other systems. Under the effect of stochastic driving, the 
system reaches a stationary state where external forces balance dissipation on average. 
These states have the invariant probability that does not respect detailed balance, and 
are characterized by non-vanishing currents of conserved quantities. In order to analyze 
spatially homogeneous stationary states, we develop a kinetic approach that generalizes 
the one known for deterministic long-range systems; we obtain a very good agreement 
between predictions from kinetic theory and extensive numerical simulations. Our approach 
may also be generalized to describe spatially inhomogeneous stationary states. We also 
report on numerical simulations exhibiting a first-order nonequilibrimn phase transition 
from homogeneous to inhomogeneous states. Close to the phase transition, the system 
shows bistable behavior between the two states, with a mean residence time that diverges 
as an exponential in the inverse of the strength of the external stochastic forces, in the limit 
of low values of such forces. 



PACS numbers: 05.20.Dd, 05.70.Ln, 05.40.-a 



1. Introduction 



Systems of particles interacting through two-body non-integrable potentials, also called 
long-range interactions, abound in nature. Common examples are plasmas interacting 
through repulsive or attractive Coulomb potential, self-gravitating systems (globular clusters, 
galaxies) involving interaction through attractive Newton potential, two-dimensional 
turbulence, and many others. In addition, there are several model systems with non- 
integrable interactions which have been studied extensively in recent years, such as spins, 
vortices in two dimensions, etc [l}|5]. 

Very often, these systems are acted upon by external stochastic forces that drive them 
out of equilibrium. Unlike systems with short-range interactions, stochastic forces in long- 
range interacting systems act coherently on all particles, and not independently on each 
particle. Consider, e.g., globular clusters being influenced by the gravitational potential of 
their galaxy, which produces a force that fluctuates along their physical trajectories. In 
addition, galaxies themselves feel the random potential of other surrounding galaxies, and 
their halos are subjected to transient and periodic perturbations, which may be due to 
the passing of dwarfs or to orbital decaying [|;6|. Dynamics of plasmas are also strongly 
influenced by fluctuating electric and magnetic fields due to the ever-changing ambiance [7| . 
In situations of stochastic driving, the systems at long times often reach a nonequilibrium 
stationary state that violates detailed balance. In such a state, the power injected by the 
external random fields balances on average the dissipation, and there is a steady flux of 
conserved quantities through the system. 

Study of nonequilibrium stationary states (NESS) is an active area of research of modern 
day statistical mechanics. One of the primary challenges in this field is to formulate a 
tractable framework to analyze nonequilibrium systems on a common footing, similar to the 
one due to Gibbs and Boltzmann that has been established for equilibrium systems [8 10 



This paper provides, to our knowledge, the first study of NESS in long-range systems with 
statistical mechanical perspectives. 

Common theoretical approaches to study isolated systems with long-range interactions 
include the kinetic theory description of relaxation towards equilibrium. In plasma physics 
and astrophysics, this approach leads to the Lenard-Balescu equation or to the approximate 



Landau equation 11 12 . One of the main theoretical results of this paper is a detailed 
development of a generalization of this kinetic theory approach to describe nonequilibrium 
stationary states in systems with long-range interactions driven by external stochastic forces, 
valid in the limit of small external stochastic fields. The nonequilibrium kinetic equation 
that we obtain describes the temporal evolution of the one-particle distribution function. In 
the limit of small external forcing, the system settles into a stationary state, in which we find 
the one-particle momentum distribution to be non-Gaussian. The predictions of our kinetic 
equation for spatially homogeneous stationary states compare very well with results of our 
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extensive A^-particle numerical simulations on a paradigmatic model of long-range interacting 
systems. Our numerical simulations also exhibit a nonequilibrium phase transition between 
homogeneous and inhomogeneous states. Close to the phase transition, we demonstrate the 
occurrence of bistability between these two types of states, with a mean residence time that 
diverges as an exponential in the inverse of the strength of the external forcing, in the limit 
of low values of such forcing. 

Similar bistable behavior has recently been observed in two-dimensional turbulence with 



stochastic forcing 13 . We believe that such phase transitions are essential phenomena for 
geophysical flows and climate, for which the two-dimensional Euler equations are a simplified 
paradigmatic model. There exists a very strong analogy between the two-dimensional Euler 
equations and the Vlasov equation relevant for leading order dynamics of the model we 



discuss in this paper 14, 15 . One of the motivations of the present work is to be able to 
study analogous phenomena in a setup for which the theory can be more easily worked out. 

In a recent letter, we reported on some of the above findings, in one of the first studies 
of non-equilibrium stationary states in systems with non-integrable potentials and driven 
by external stochastic fields (16]. The aim of this paper is to present a detailed derivation 



of the results given in Ref. 16 , as well as to report on additional empirical results, more 



specifically non-equilibrium phase transitions. 

We note that Ref. 17 presents a computation of the same kinetic equation as the one 



we describe in this work and in 16 . Nevertheless, the result is different. 

The structure of the paper is as follows. In the following section, we define the dynamics 
we are going to consider of a long-range interacting system driven by external stochastic 
forces. We also discuss the paradigmatic example of the Hamiltonian mean-field (HMF) 
model. In section [3| we discuss the methods we adopt to analyze the dynamics. In particular, 
we give a detailed derivation of the kinetic theory to study spatially homogeneous stationary 
states of the dynamics. We describe the numerical simulation scheme that we employ to study 
the dynamics, specifically, to check the predictions of our kinetic theory. This is followed by 
a discussion in section |4] of the results obtained from the kinetic theory, and their comparison 
with numerical simulation results for spatially homogeneous stationary states. In section [5} 
we discuss the results of numerical simulations of spatially inhomogeneous stationary states. 
We report on the very interesting bistable behavior in which the system in the course of 
its temporal evolution switches back and forth between homogeneous and inhomogeneous 
states, with a mean residence time that we show to be diverging as an exponential in the 
inverse of the strength of the external stochastic forcing, in the limit of low values of such 
forcing. We close the paper with concluding remarks. Some of the technical details of our 
computation are collected in the four appendices. 
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2. Long-range interacting systems driven by stochastic fields 



2.1. The model 

Consider a system of particles interacting through a long-range pair potential, and 
described by the Hamiltonian 

N 2 1 ^ 
i=l i,j=l 

Here, qt and pi are, respectively, the coordinate and the momentum of the i-th particle, and 
v{q) is the two-body interaction potential. We take the particles to be of unit mass. In this 
paper, for simplicity, we regard g^'s as scalar periodic variables of period 27r; generalization 
to qi G M", with = 1, 2 or 3, is straightforward. 

In plasma physics, the typical number of particles interacting with one particle is given 
by the coupling parameter T = nXjj, where n is the number density, and is the Debye 
length. It is then usual to rescale time such that the inverse of F multiplies the interaction 
term jll]. In self-gravitating systems, the dynamics is dominated by collective effects, so that 
it is natural to rescale time in such a way that the parameter multiplies the interaction 
potential 18 . These facts explain the rescaling of the potential energy by in Eq. 



called the Kac scaling in systems with long-range interactions 19 . We emphasize that no 
generality is lost in adopting the Kac prescription. 

We perturb the system ([T]) by a statistically homogeneous Gaussian stochastic field 
F{q, t) with zero mean, and variance given by 

{F{q,t)F{q',t'))=C{\q-q'\)5{t-t'). (2) 

The resulting equations of motion for the i-th particle are 

dH ^ . dH ^r., ^ 

qi = — , and = api + y/a F{qi, t). (3) 

opi oqi 

The property that the Gaussian fields F{qi,t) are statistically homogeneous, i.e., the 

correlation function C depends solely on \q — q'\, is consistent with any perturbation that 

respects space homogeneity. Such a property is necessary for the discussions later in the 

paper on the kinetic theory approach to describe spatially homogeneous stationary states 

of the dynamics Note that C{q) is the correlation, so that it is a positive-definite 

function [20j, and its Fourier components are positive: 

Cfc = — / dq Ciq)e-''"' > 0; c.^ = c^, C{q) = Cq + 2 V cos(A:g). (4) 

We find it convenient to use the equivalent Fourier representation of the Gaussian field F{q, t) 
as follows: 

oo 

F{q, t) = ^0 Mt) + [cos(A;g) X,,(t) + sin(A;g) Yk{t)] , (5) 

k=l 



where Xo(t), Xk{t) and Yk{t) are independent scalar Gaussian white noises satisfying 

{Xk{t)Xkit')) = 6k,k'S{t-t'), (6) 

{Yk{t)Ykit')) = 6k,k,6{t-t'), (7) 

{Xk{t)Yk'{t')) =0. (8) 

For stochastic dynamics of the type Eq. (|3|, general mathematical results allow to prove 
that the dynamics is ergodic (see, for instance, [2l]). 



Using the Ito formula 22 to compute the time derivative of the energy density e = H/N, 



and averaging over noise realizations give 

(I) ^ <2-> = (^) 

where k = J^iLiPl / i'^^) is the kinetic energy density. On integration, we get for 
homogeneous states for which e = n that 

(kit)) = ((MO)) - ^) e-- + (10) 

The average kinetic energy density in the stationary state is thus (k)^^ = C(0)/4. We define 
the kinetic temperature of the system as 

{'^)ss = (11) 

as a result, we have 

r = ^M. (12) 

2 ^ ' 

Let us note that in the dynamics (|3]), fluctuations of intensive observables due to 
stochastic forcing are of order ^/a, while those due to finite-size effects are of order 1 / 
Moreover, the typical timescale associated with the effect of stochastic forces is l/a, as 



is evident from Eq. (10), while the one associated with relaxation to equilibrium due to 
finite-size effects is of order A^, see jl]|2|. 

Our theoretical analysis to study the dynamics ([s]) by means of kinetic theory is valid 
for any general two-particle interaction potential v{q). However, in order to perform simple 
numerical simulations with which we may check the predictions of the kinetic theory, we 
specifically make the choice v{q) = 1 — cosg, that defines the stochastically-forced attractive 
Hamiltonian mean-field (HMF) model, as detailed below. 

2.2. A specific example: The stochastically-forced HMF model 

The Hamiltonian mean-field (HMF) model is a paradigmatic model to study long-range 
interacting systems. The model describes particles moving on a circle under deterministic 
Hamiltonian dynamics, and interacting through the interparticle potential v{q) = 1 — cosg 



23 ,24 . It displays many features of generic long-range interacting systems, e.g., existence of 
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quasistationary states [T,24 . In equilibrium, the model has a second-order phase transition 
from a high-energy spatially homogeneous phase to a low-energy inhomogeneous phase at the 
energy density Cc = 3/4, corresponding to the critical temperature = 1/2. In a system of 
particles, the degree of spatial inhomogeneity at time t is measured by the magnetization 
variable m{t), defined as 



m{t) 



1 

N 



\ 



N 



cos qi + 



TV 

^ sin Qi 



(13) 



In the thermodynamic limit — )■ oo, the magnetization in the steady state decreases 
continuously as a function of energy from one to zero at the transition energy Cc, and remains 
zero at all higher energies. When forced by the stochastic forces F{qi,t) resulting in the 
dynamics ([3]), we call the corresponding model the stochastically-forced HMF model. We 
note for later purpose that the Fourier transform of the HMF interparticle potential is, for 
k 0, Vk = — [Sk,i + Sk,-i] /2, where 6k,i is the Kronecker delta. 



3. Methods of analysis 

3.1. Kinetic theory for homogeneous stationary states 

Here, we develop a suitable kinetic theory description to study the dynamics ^ in the 
joint limit — )■ oo and a — )■ 0. While the first limit is physically motivated on grounds 
that most long-range systems indeed contain a large number of particles, the second one 
allows us to study stationary states for small external forcing. Moreover, for small a, we will 
be able to develop a complete kinetic theory for the dynamics. For simplicity, we discuss 
here the continuum limit Na ^ 1, when stochastic effects are predominant with respect to 
finite-size effects. The generalization of the following discussion to the cases Na of order 
one and Na ^ 1 is straightforward, as pointed out at the end of this subsection. For the 
development of the kinetic theory, we assume the system to be spatially homogeneous; a 
possible generalization to the non-homogeneous case will be discussed in the conclusions of 
the paper. 

As a starting point to develop the theory, we consider the Fokker-Planck equation 
associated with the equations of motion ([3]). This equation describes the evolution of the 
A^-particle distribution function fN{qi,---,qN,Pi,---,PN,t), which is the probability density 
(after averaging over noise realizations) to observe the system with coordinates and momenta 
around the values {qi,Pi}i<i<N at time t. This equation can be derived by standard 
methods [22]; we have 



dt 



dfw , d{apifN] 
-Pi— \- 



dqi 



dpi 



1 

2N 



dqi 



d_ 
dpi 



d_ 

dp J J 



f. 



N 



6 



N 



dpidpj 



(14) 



In Appendix A , by applying the so-called potential conditions 25 for the above Fokker- 
Planck equation, we prove that a necessary and sufficient condition for the stochastic 
process ^ to verify detailed balance is that the Gaussian noise is white in space, that 
is, Ck = c for all k. This condition is not satisfied for a generic correlation function C{q), in 
which case, the steady states of the dynamics are true nonequilibrium ones, characterized by 
non-vanishing probability currents in configuration space, and a balance between external 
forces and dissipation. 

Similar to the Liouville equation for Hamiltonian systems, the iV-particle Fokker-Planck 
equation ( 14 ) is a very detailed description of the system. Using kinetic theory, we want to 
describe the evolution of the one-particle distribution function 



N 



f{zi,t) 



YldZifN{zi,...,ZN,t), (15) 
^ i=2 

where we have used the notation zi = {qi,Pi). We note that with this definition, the 
normalization is J dzf{z,t) = 1. Substituting in the Fokker-Planck equation ( 14 ) the reduced 
distribution functions 

N 



fs{Zi, ...,Zs,t) 



{N-s)\N' 



Zn, t), 



(16) 



i=s+l 



and using standard techniques [26], we get a hierarchy of equations, similar to those of the 
Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) hierarchy, as follows: 

^ dfs 1 dv{qi - qj) dfs ^ d 

dt 



dq, N 



d'fs 



dpidpj 



dqi 

s 



dpi 



dzs+i—v{qi - g^+i)— — 
dqi dpi 



(17) 



for s 



i,j=l ' i=l 

N — 1. In this paper, we use both the notations dh/dq and h' to denote the 



derivative of a function h. With a slight abuse of the standard vocabulary, we will refer to 
Eq. (17) as the BBGKY hierarchy equation. 

Now, as is usual in kinetic theory, we spht the reduced distribution functions into 
connected and non-connected parts, e.g., 

f2iZi,Z2,t) = f{zi,t)f{z2,t) +g{zi,Z2,t), (18) 

and similarly, for other /^'s with s > 2. In Appendix B we show that the connected part 
g{zi, Z2, t) of the two-particle correlation is of order a, so that we may write 

f2izi,Z2,t) = fizi,t)f{z2,t) + ag{zi,Z2,t), (19) 
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where g is of order unity; more generally, the connected part of the fc-particle correlation 
is of higher order, with respect to a, in the small parameters a and Then, to close 

the BBGKY hierarchy, we neglect the effect of the connected part of the three-particle 
correlation on the evolution of the two-particle correlation function. This scheme is justified 
at leading order in the small parameter a, and is the simplest self-consistent closure scheme 
for the hierarchy while taking into account the effects of the stochastic forcing. With our 
assumption that the system is homogeneous, i.e., / depends on p, and g depends on \qi — q2\, 
Pi and p2 only, the first two equations of the hierarchy are then 

d 



df d 



'^dp ' ^^^^P^'"'^^'^^3{qi,p,Put), 



and 



+ Lyg + L) 'g 



dt 



C(|gi-g2|)/'(Pi,t)/'(P2,t), 



(20) 



(21) 



where L^P and Ly are the Vlasov operators linearized about the one-particle distribution /, 



'(2) 



and acting, respectively, on the first pair (gi, pi) and on the second pair (g2, P2) of the function 
g = g{qi,Pi,q2,P2,t). Exphdtly, for a function h of {q,p,t), the expression for the linear 
Vlasov operator Lfh is 

dh 



Lfh{q,p,t) = p— - f'{p,t) J dqidpiv'iq- qi 



)h{qi,pi,t), 



22 



so that we have for L^^^ 



f 3. 

Lfg{qi - q2,Pi,P2,t) 



Pi 



dg 
dqi 



f'iPi,i) j dqsdpsv'iqi - q3)g{q3 - q2,P3,P2,t). 

(23) 



Ljr g is obtained from Eq. (23) by exchanging the subscripts 1 and 2. 

To obtain from Eqs. (20) and (21) a single kinetic equation for the distribution 
function /, we have to solve Eq. (21) for g as a. function of / and plug the result into 
the right hand side of Eq. (20). Because the two equations are coupled, this program is 
not achievable without making further simplifying assumption. Nevertheless, we readily see 
from these equations that the two-particle correlation g evolves over a timescale of order 
one, whereas the one-particle distribution function f{p,t) evolves over a timescale of order 
1/a. We may then use this timescale separation and compute the long-time limit of g from 
Eq. (21) by assuming / to be steady in time; this is the equivalent of the Bogoliubov's 
hypothesis in the kinetic theory for isolated systems with long-range interactions. Note 
that for this timescale separation to be valid, we must also suppose that the one-particle 
distribution function f{p,t) is a stable solution of the Vlasov equation at all times. Indeed, 
if this is not the case, it can be shown that g diverges in the limit t — ?• 00 27 . The physical 
content of this hypothesis is that the system slowly evolves from the initial condition through 
a sequence of quasistationary states to the final stationary state. 
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Because we assume the system to be homogeneous in space, it is useful to Fourier 



transform Eqs. (20) and (21) with respect to the spatial variable; we get 



df d 



d f 

-2nia ^ ^^'^^ J ^P' 9k{p,p' ,t), 



and 



dgk 
dt 



+ Lflgk + Lfl^Qk^ {pi,P2,t) = Ckf\pi)f\p2), 



(24) 



(25) 



where gk{pi,P2, t) is the Fourier transform of g{q, Pi,P2,t) with respect to the spatial variable, 
and Vk is the k-th Fourier coefficient of the pair potential v{q). The explicit expression for 
the k-th Fourier component of the linear Vlasov operator Lf ^ acting on a function h{p) is 



Lf^kh){k,p) = ikph{p) -2TTikvkf'{p) / dp h{p'). 



(26) 



One has analogous expressions for L^^]. and for Ly_j^. We readily see that L^^, = L 



(2) 



f-k- 



From the right hand side of Eq. (24), we see that to obtain a single kinetic equation. 



we need only the Fourier transform gk{p,p',t), more specifically, its integral with respect to 
the second momentum variable p'. Actually, it can be shown that gk{p,p', t) does not have a 
well-defined time-asymptotic (it converges only in the sense of distribution), while its integral 



with respect to p' does have; this is connected to the mechanism of Landau damping 27 



The structure of Eqs. (20) and (21), or, equivalently, of Eqs. (24) and (25) is very familiar 



in kinetic theories; we refer the reader to 27 for a general discussion. Equation (21) 



or, equivalently, Eq. (25), is called the Lyapunov equation for the two-point correlation 
of a stochastic variable described by an Ornstein-Uhlenbeck process. However, there 
is a difference from the standard finite-dimensional case 22 in that in our case, Lf is 



a linear infinite-dimensional operator acting on a functional space, instead of being a 
finite-dimensional one, i.e., a matrix. This makes it non-trivial to compute the long-time 



asymptotic of the right hand side of Eq. (20 ), where g is the solution of Eq. (21 ) with / steady 



in time. A possible way to achieve this goal is to follow the derivation of the Lenard-Balescu 



equation from the BBGKY hierarchy, as may be found in the Appendix A of Ref. 11 . For 



explicit technical details, see 27 , in which the method to solve the Lyapunov equation in a 



general manner is discussed, and, subsequently, applied to the derivation of kinetic theories 
for long-range interacting systems and two-dimensional turbulence models. 

In the present case, the linear transform of the stationary solution of the Lyapunov 



equation, Eq. (25), which is needed to compute the right hand side of Eq. (24), can be 



written (see Ref. [27)) in the frequency space as 



dpig^[f]iP,Pi) = hra / dpigkip,pi,t) (27) 
= lj^d^ {RfA^)b)ip) J dp' (i?;,_,(-c^)6*)(p'), (28) 



where F is a contour which passes above all the poles of (^Rf^k{(^)bj , and Rf^k{^) is the 
resolvent operator, defined as 

Rfj^{uj) = {-iuj + Lf^kY\ (29) 

while b(j),t) = y/ckf'{p,t). The discussion of the explicit form of the resolvent operator is 
a standard topic in plasma theory, and involves the phenomenon of Landau damping; we 
refer the reader to classical references for this result, for example, [Tl|[l2||28]. Its action on 
a function h, defined for u such that Im(ci;) > 0, is 



Rf^k{(^)h){p) 



■ 2Tiikvk . 

h{P) + /, J [P) 



dp' 



Kp') 



—iu + ikp 

where e{k,u) is the dielectric function, which for Im(ci;) > is given by 

f'ip) 



—iu + ikp' 



e{k, uj) 



27civkk / dp 



(30) 



(31) 



-ico + ikp 

and by its analytic continuation for u when lm{u)) < 0. Both the resolvent operator and the 
dielectric function are defined for G M by their analytic continuation, which will still be 
denoted by the same symbols. 

Now, inserting Eq. (30) into Eq. (27), and with some calculations whose details will be 



(32) 



reported in |27|, we get the kinetic equation 



where 



df Jjpf) ^ d 
dt dp dp 



m{p) = ^cio) 



D[f] 



dp 



VkCk 



k=l 



dpi 



+ 



\e{k,kp)\'^ \e{k,kpi)\'^ \ Pi — p 



-/'(Pi,t)(33) 



We recall that Vk is the k-th Fourier coefficient of the pair potential v{q), the quantity Ck 
is defined in Eq. (ji), e is the dielectric function defined in Eq. (31), and J* indicates the 
Cauchy integral or Principal Value. 

The kinetic equation (32) is the central result of the kinetic theory developed in this 
paper. It has the form of a non-linear Fokker-Planck equation 25 because the diffusion 
coefficient D[f]{p) is itself a functional of the one-particle distribution function /. The 
linear part of the diffusion coefficient (1/2)C(0) is the mean- field effect of the stochastic 
forces, whereas the effect of two-particle correlation is encoded in the non-linear part. In 
the next section, we describe how we use this kinetic equation to get information about the 
nonequilibrium stationary states of the dynamics. 

In the foregoing, we discussed the kinetic theory in the limit Na ^ 1. The extension to 
the general case is straightforward: Because of the linearity of the equations of the hierarchy 
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(20) and (21 ), the finite- and stochastic effects give independent contributions. The kinetic 
equation at leading order of both stochastic and finite-size effects is 

j^ = QM] + QN[f], (34) 



where Qa is the operator described in Eq. (32), and Qn (of order 1/A^) is the Lenard- 



Balescu operator. For instance, in the case Na ^ 1 and in dimensions greater than one, 
the operator Qn is responsible for the relaxation to Boltzmann equilibrium after a timescale 
of order A^, whereas the smaller effect of Qa selects the actual temperature after a longer 
timescale of order 1/a. 

3.2. Numerical simulations 



1 




200 400 600 800 1000 



Time t 



Figure 1. Magnetization as a function of time, obtained from numerical simulation of the 
stochastically-forced HMF model with N ~ 1000, a = 0.01 at kinetic temperature T = 0.25, 
and with ci — C2 — C'l, ~ . . . ~ ciq ~ C(0)/20, Cfc>ii = 0, where C(0) — 2T. The values of 
the integration step size At used are marked in the figure. The data are obtained by using 
the integration algorithm described in section |3.2[ That the magnetization plots collapse 
onto one curve shows the stability of our algorithm with respect to variation in At. We 
have checked that the final value of the magnetization matches with the prediction from 
equilibrium statistical mechanics. 



Here we describe how we may simulate the dynamics ^ by means of a numerical 
integration scheme. To simulate the dynamics over a given time interval [0 : 7~], choose 
a time step size At, and set t„ = nAt as the n-th time step of the dynamics. Here, 
n = 0, 1, 2, . . . , Nt, where Nt = T/At. In our numerical scheme, at every time step, we first 
discard the effect of the noise and employ a fourth-order symplectic algorithm to integrate 
the deterministic Hamiltonian part of the dynamics 29 . Subsequently, we add the effect of 
noise and implement an Euler-like first-order algorithm to update the dynamical variables 
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H} Specifically, one step of the scheme from t„ to tn+i = t„ + At involves the following 
updates of the dynamical variables for i = 1,2, . . . , N: For the symplectic part, we have, for 



m 



Pi[tn+ 



mAt 
mAt 



Pi(tn + 7^ — ) + b{m)At 



OH,, ( (m-l)At\,/ 



dq, 



( m/\t\ ( 



im - l)At 



+ a{m)At Pi {tn + 



mAt 



where the constants a(m)'s and b{mys are given in Ref. [29]. At the end of the update (35), 
we have the set {qi(tn+i) , Piitn+i)} ■ Next, one includes the effect of the stochastic noise by 
leaving gi(tn+i)'s unchanged, but by updating pj(t„+i)'s as 



1 - aAt 



Piitn+l ) Piitn+l] 
Nr 

+ J2 V^k{AX^''\tn+i) COS (kq,{tn+i)) + Ar(^-)(t„+i) sin (kqi{tn+l] 



k=l 



.(36) 



(35) 



Here AX^^^ and AF^^') are Gaussian distributed random numbers with zero mean and unit 
variance. The outcome of implementing this mixed scheme for the stochastically-forced HMF 
model is shown in Fig. [T| where one may observe consistent results with respect to change 
of At over a wide range of values. In numerical simulations reported later in the paper, we 
exclusively used this mixed scheme to simulate the dynamics (|3|). 



4. Predictions of the kinetic theory and comparison with simulations 



We now focus on how to obtain from the kinetic equation (32) predictions for the 



nonequilibrium stationary states of the system. According to Eq. (32), 1/a is only a 
timescale; thus, at leading order in a and except for a time rescaling, the parameter a 
does not affect the time evolution of the system. This statement holds also beyond the 
leading order for what concerns the evolution of the kinetic energy; its evolution may be 
obtained directly from the equations of motion as discussed in section |2| and can also 



be obtained from the kinetic equation (32), as detailed in Appendix C For the evolution of 



other observables, there will be corrections at higher orders in a. 

As previously discussed, the dynamics of the system does not respect detailed balance 
if the forcing is not white in space. At the level of the kinetic equation, by inspecting the 



definition of the diffusion coefficient, Eq. (33 ), we see that the effect of correlations induced by 



the stochastic forces is modulated by the Fourier component Vk of the interparticle potential. 

I We found that an Euler-like first-order scheme afone is unstable with respect to not-too-small At, in the 
sense that one obtains different magnetization profiles as a function of time t = t„At. The situation gets 
worse for small a, when one needs to use very small At to obtain consistent results. Therefore, for faster 
and efficient simulation, we adopted the "mixed" scheme described in the text. 
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Figure 2. (a) Kinetic energy density (k), and (b) (p^) as a function of at, at kinetic 
temperature T = 0.75, with modes 1—2 excited with amphtudes satisfying ci — C2 — C(0)/4, 
where C(0) = 2T. The data for different and a values are obtained from numerical 
simulations of the stochastically- forced HMF model with At = 0.01, and involve averaging 
over 50 histories for N — IC* and lO'^ histories for N = 10'^. The data collapse implies that 
a is the timescale of relaxation to the stationary state. The inset shows the data without 



time rescaling by a. Similar plots for different parameter values were reported in Ref. 16 
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Figure 3. Stationary momentum distribution f{p), on (a) linear, and (b) semi-log scales, 
for a = 0.001, and 0.01 at kinetic temperature T — 0.75. The plots correspond to modes 
1 and 2 excited with amplitudes satisfying ci — C2 — C(0)/4, where C(0) = 2T. The 
data denoted by crosses and squares are results of iV-body simulations of the stochastically- 
forced HMF model with N = 10000, At = 0.01 and 1000 independent reahzations of the 
dynamics, while the red continuous lines refer to the theoretical prediction from the kinetic 
theory. For comparison, the black broken line shows the Gaussian distribution with the 
same kinetic energy (stationary state of the stochastically-forced HMF model at T = 0.75, 
Co = ci = 0, C2 = 0.75, Ck>3 = 0). 



Then, taking the forcing spectral amphtudes Ck different from zero if and only if Vk = 0, 
the non-linear part of the diffusion coefficient vanishes. On the other hand, taking Ck ^ 
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Figure 4. Diffusion coefficient D[f]{p) for the stationary momentum distribution f{p) at 
kinetic temperature T = 0.75, with cq = Cfc>3 0, and either (i) ci ^ C2 — 0.375, or, (ii) 
ci = 0,C2 = 0.75. 



for the modes for which Vk ^ leads to a diffusion constant which has a non-vanishing 
non-hnear part. To be concrete, let us discuss these two scenarios in the context of the 
stochastically-forced HMF model. 

Since the Fourier transform of the HMF interparticle potential is, for k ^ 0, Vk = 
— [Sk,i + ^k-i] /2, it follows that only the stochastic force mode with wave number k = 1 
contributes to the non-linear part of the diffusion coefficient; all the other stochastic force 
modes result in only a mean- field contribution through the term C(0). Thus, for the case 
Ci 7^ 0, the two relevant parameters that dictate the evolution of the stochastically-forced 
HMF model by the kinetic equation (32) with a non- linear diffusion coefficient are C(0) 
and Ci. From Eq. (12), since C(0) is related to the kinetic temperature T, we take T and 
Ci to be the two relevant parameters. From Eq. (11), we know that 2T equals the kinetic 
energy in the final stationary state. Also, Eq. ^ implies that ci < C(0)/2. 

If however ci = 0, then, at leading order in a, the dynamics of the system is described 
by a linear Fokker-Planck equation; this equation is the same as the one which describes the 
HMF system when coupled to a Langevin thermostat, studied in 30,31 . This means that for 
this particular choice of the parameters, the detailed balance is broken for the dynamics, but 
this feature cannot be seen in the kinetic theory, being an effect at a higher order in a. In this 
case, the homogeneous stationary states of the kinetic equation have Gaussian momentum 
distribution f{p). As has been studied thoroughly in the context of canonical equilibrium 
of the HMF model, these states are stable for kinetic energies greater than 1/4, i.e., for 
C{0) > 1. 

Except for the special case of ci = 0, the stationary velocity distribution of the kinetic 
equation (32) is in general not Gaussian. This can be seen semi-analytically by observing 
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that the Gaussian distribution function 

= A exp(-/3p2) , A=\^, /3 = -^, (37) 

V vr Zi 

with /3 chosen such that the value of the kinetic energy is the one selected by T, solves the 
linear Fokker-Planck equation with the diffusion coefficient given by 

Dmf = T. (38) 



To prove that the Gaussian distribution function is not a stationary solution of Eq. (|32j), 
we have to prove that the contribution to df/dt from the non-linear part of D[f]{p) in 
Eq. (33) does not vanish. This result can be proven with an asymptotic expansion [12 
for large momenta of the integrals which appear in the diffusion coefficient. We report the 
straightforward computation in Appendix D From the same analysis, one can deduce that, 
even though the distribution function is not Gaussian, its tails are Gaussian. 

On the basis of the above discussions, we expect that for values of T and Ci such that 
T > 0.5 and ci <^ 2T, the stationary states will be close to homogeneous states with Gaussian 
momentum. In order to locate the actual stationary states of the kinetic equation, we have 
devised a simple numerical scheme, based on the observation that a linear Fokker-Planck 
equation whose diffusion coefficient D{p) is strictly positive admits a unique stationary state 



fss{p) = Aexp 



dp ^ 



(39) 



D{p') 

For a given distribution fn{p), we compute the diffusion coefficient Dn{p) through Eq. (|33 



and then fn+i using D„ and Eq. (39). This procedure defines an iterative scheme. Whenever 
convergent, this scheme leads to a stationary state of Eq. (32). Each iteration involves 
integrations, so that we expect the method to be robust enough when starting not too 
far from an actual stationary state. However, we have no detailed mathematical analysis 
yet. Implementing this iterative scheme, we observed that the distribution to which the 
scheme converges is independent of the initial distribution /q. Moreover, the convergence 
time is exponential in the number of steps n whenever T is not too close to loss of stability 
of /oo with respect to the linear Vlasov dynamics; in practice, we are able to get reliable 
results for T > 0.65. 

In order to check the theoretical predictions discussed above, we performed numerical 
simulations of the stochastically forced HMF model. Figure [2] shows the evolution of 
the kinetic energy and (p^) = (1/A^) XliliPf) where they have been compared with our 
theoretical predictions. In the case of (p'*), we have compared the long-time asymptotic 
value with the kinetic theory prediction for the stationary state, computed numerically by 
using the iterative solution for the stationary distribution. The figure illustrates a very good 
agreement between the theory and simulations. For a more accurate check of the agreement, 
we have obtained the stationary momentum distribution from both A^-body simulations and 
the numerical iterative scheme discussed above. A comparison between the two, shown in 
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Fig. [3} both on linear and semi-log scales, shows a very good agreement between theory and 
simulations. In this figure, we also show the Gaussian distribution with the same kinetic 
energy, to illustrate the point that the stationary momentum distribution of the system is 
far from being Gaussian. The diffusion coefficient D[f](p) is shown in figure |4j 

In passing, let us remark that, with an iterative scheme analogous to the one described 
above, one could have also obtained the full time evolution f{p,t) that obeys the kinetic 



equation (32). However, we will not address this point here. 

We also note that while a linear Fokker-Planck equation with non-degenerate diffusion 
coefficient can be proven to converge to a unique stationary distribution |25|, this is not 



true in general for non-linear Fokker-Planck equations like Eq. (32). We expect that if 
the dynamics is not too far from detailed balance, the kinetic equation will have a unique 
stationary state. Far from equilibrium, the kinetic equation could lead to very interesting 
dynamical phenomena, like bistability, limit cycle or more complex behaviors. The main 
issue is then the analysis of the evolution of the kinetic equation. Although some methods to 



study this type of equation exist 32 , we have only the numerical iterative scheme described 
above to provide some preliminary answers. A more rigorous mathematical analysis is left 
for future studies. 

5. Nonequilibrium phase transition and collapse 

Until now, we have considered homogeneous stationary states of the dynamics ([3]), and 
have discussed a kinetic theory to analyze them. Although our theory can in principle be 
extended to include inhomogeneous stationary states, its actual implementation to get, e.g., 
the single-particle distribution, would require more involved computations than the one we 
encountered for homogeneous states. In order to get preliminary answers, we have performed 
extensive numerical simulations of the dynamics in the context of the stochastically-forced 
HMF model. Our specific interest is to know about how the magnetization behaves as the 
kinetic temperature is reduced from high values. 

In the case when the stochastic forcing respects detailed balance (i.e., when the noise 
spectrum is fiat and all modes are excited), the stochastically-forced HMF model reduces 



to the Brownian mean-field (BMF) model studied previously 31 . Here, we know that the 
system settles into an equilibrium state in which it exhibits a second-order phase transition 
at the kinetic temperature T = T^. = 1/2: on increasing T from low values, the magnetization 
decreases continuously to zero at Tc and remains zero at higher temperatures. In the 
following, we excite only a limited number of modes Nji, but the amplitudes of all excited 
modes are equal (c^ equals c for all k < N^, and is zero otherwise, where the constant 
c is related to the temperature). Figure |5|^a) shows that with Nr = 50, one reproduces 
very well the equilibrium profile of the magnetization as a function of temperature. On 
reducing the value of N^, the system is driven more and more out of equilibrium. Indeed, 
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Figure 5. Numerical simulation results for magnetization in the stochastically-forced 
HMF model as a function of adiabatically-tuned kinetic temperature T; the different plots 
correspond to different number of modes excited in the spectrum with amplitudes satisfying 
C(0) = Co + 2X;^=i Cfc, where C(0) = 2T, the index k = 1,2,...,Nr denotes the mode 
number, while Nfi is the total number of excited modes with fc 7^ 0. In all cases, the modes 
excited were chosen to have equal amplitudes, with cq = 0, N = 5000, a = 0.01, At = 0.01, 
while the tuning rate for T is 10~^. It may be noted that forcing equally a large number 
of modes {'^ 50) reproduces the equilibrium magnetization profile as illustrated by the 
match with the analytical equilibrium solution in the panel (a). In panel (b), the first-order 
nonequilibrium phase transition is marked by the vertical dashed line. In panel (c), besides 
the first-order transition, we also show the dynamical transition to the collapsed state by 
the vertical dashed dotted line. In panel (d), the nonequilibrium phase transition and the 
dynamical transition almost coincide, and we just show the latter one by the vertical dashed 
dotted line. 



Fig. [sj^b) shows that with Nr = 7, the magnetization profile changes; in particular, it 
develops a discontinuity around a temperature Ttrans ~ 0.49, reminiscent of a first-order 
phase transition. The transition temperature is denoted by the vertical dashed line. With 
Nji = 3, Fig. [sj^c) shows that the discontinuity gets more pronounced, and Ttrans is now 
shifted to a higher value (denoted again by the vertical dashed line). A new feature appears 
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in this plot, namely, at a temperature T^yn ~ 0.4, the magnetization attains the maximal 
value of unity, which it retains for all lower temperatures. This value of unity corresponds 
to a state in which the particles are very close to one another on the circle, thus defining 
a "collapsed" state. We found that this state, as well as the transition to it, persist on 
changing the system size A^. 

Now, it is known that trajectories of ensembles of dissipative dynamical systems forced 
by the same realization of a stochastic noise converge to a single one 33 34 . These attracting 
trajectories are referred to as the ones due to the so-called stochastic attractor. Although 
we did not perform a detailed characterization of the collapse in our model, we believe that 
the phenomenon is related to stochastic attractors. 

Coming back to Fig. [sj^c), we see that for temperatures T^yn < T < Ttrans, the 
magnetization shows strong fluctuations. Reducing the number of excited modes to a single 
one, namely, to the one that coincides with the Fourier mode of the HMF potential, it seems 
from Fig. [sj^d) that only the dynamical transition to the collapsed state at a temperature 
Tdyn ~ 0.66 persists. 

The hint that the nature of the phase transition at Ttrans is of first-order comes from 
the hysteresis plots of Fig. |6| To obtain these plots, one monitors the magnetization while 
tuning adiabatically the kinetic temperature across Ttrans from higher to lower values and 
back to complete a full cycle. As is evident from Fig. [6} the observed hysteresis is between 
the collapsed state and the zero-magnetization state. In principle, it should be possible to 
observe a hysteretical behavior between the magnetized and the zero-magnetization state. To 
achieve this in simulations involving adiabatic tuning of temperature, one should not allow 
the system to make the transition to the collapsed state, which requires conditions close to 
those that ensure detailed balance. However, a possible drawback of this method is that 
closeness to detailed balance might lead to narrow hysteresis loops. Moreover, the adiabatic 
tuning of temperature should not be very slow, as otherwise one observes bistability instead 
of the hysteresis. All these factors make the observation of hysteresis between the magnetized 
and the zero-magnetization state difficult to observe numerically; further explorations of this 
will be the subject of future investigations. 

In order to explore further the region in Fig. [sj^c) close to Ttrans, and to ascertain the 
nature of the phase transition at Ttrans, we fix the value of the temperature to be T = 0.53, 
and monitor the magnetization as a function of time. The time series of the magnetization 
is shown in Fig. [7]^a), in which one observes clear signatures of bistability, whereby the 
system switches back and forth between homogeneous (m ^ 0) and inhomogeneous (m > 0) 
states. In addition, we show in Fig. [7](b) the distribution of the magnetization around 
the phase transition temperature: the distribution is bimodal with a peak around a zero 
value and another one around a positive value. When decreasing the temperature across 
the phase transition region, we clearly see that the peak heights of the distributions of the 
magnetization at the zero and non-zero values interchange. These two features, together 
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Figure 6. Numerical simulation results for magnetization in the stochastically-forced HMF 
model as a function of adiabatically-tuned kinetic temperature T for two different values of a. 
In each case, the modes 1 — 3 are excited with amplitudes satisfying ci = C2 = C3 = C(0)/6, 
where C(0) = 2T . In all cases, N = 5000, At = 0.01, while the tuning rate for T is lO^^. 
The grey points correspond to the case when the temperature is decreased from high values, 
while the black points correspond to the case when the temperature is increased from low 
values. 



with the hysteresis plots of Fig. [6| support the first-order nature of the transition around 
T^trans which cau be estimated from Fig. [7|^b) to be Titans ~ 0.532. 

From Fig. [7| it is clear that the system has two well separated attractors, corresponding 
to the homogeneous and inhomogeneous states. A question of immediate interest is: How 
long does the system stay in one state before switching to the other? Let us define the 
residence time as the time the system stays in one state before it switches to the other. 
In the limit of low noise level a, there is a clear separation between the natural dynamical 
time and the typical residence time, as is evident from Fig. [Tf^a). As a result, one may 
conjecture that two successive switching events are statistically independent of one another. 
In case such a conjecture holds for our model, the residence time statistics will be a Poisson 
process, characterized solely by the probability per unit time, A+, of switching from the 
inhomogeneous state to the homogeneous state, and the probability per unit time, A_, for 
the reverse switch. The distribution of residence time r in each phase is then exponential: 

P±(r) = -^exp(-A±r), (40) 

so that the average residence times in the two states are r^^* = Such an exponential 
form of the residence time distribution is verified from our simulation data displayed in Fig. 
[8| Note that generating such a plot requires running simulations of the dynamics for long 
enough times so that the magnetization switches back and forth between the two states a 
sufficient number of times, and one has good statistics for the residence times. For low values 
of a, such as those used in Fig. [8| this was often not feasible due to very long simulation 
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Figure 7. (a) Numerical simulation results for magnetization in the stochastically- 
forced HMF model as a function of time at kinetic temperature T = 0.53, with N = 
5000, a — 0.005, Ai = 0.01, and with modes 1 — 3 excited, whose amplitudes satisfy 
ci = C2 — c-3 ~ C(0)/6, where C(0) = 2T. The figure shows clear signatures of bistability in 
which the system during the course of evolution switches back and forth between spatially 
homogeneous (m ~ 0(0)) and inhomogeneous (m ~ 0(1)) states, (b) Distribution Prob{m) 
of the magnetization m as a function of T at a fixed value of a = 0.01. The data are obtained 
from numerical simulation results similar to (a) for magnetization in the stochastically-forced 
HMF model, with N = 5000, Ai — 0.01, and with modes 1 — 3 excited, whose amplitudes 
satisfy ci = C2 = C3 = C(0)/6, where C(0) = 2T. 



Modes 1-3 





-14 






Ot=0 .00625 • 




-15 






ft=0 .00875 X 




-16 








+ 
CL 


-17 
-13 


. y. * 

• 

X 


• 


• 




-19 






• 




-20 






• 




-21 









150000 300000 

T 

Figure 8. Distribution of the residence time r in the inhomogeneous state, for two values 
of a. The data are obtained from simulations with modes 1 — 3 excited, whose amplitudes 
satisfy ci = C2 = C3 = C(0)/6, where C(0) = 2T. Here, the kinetic temperature T = 0.53, 
while At = 10^2 iV = 5000. 



times. This results in bad statistics, and hence, the form of the plot displayed in Fig. [8} 
which, though good, may be improved upon by running longer simulations. We conclude 
that our conjecture of two successive jumps being independent holds for our model, and that 
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the average residence time fully characterizes the switching process for small enough a. 

We now discuss how the residence times depend on the system parameters, in particular, 
on a. For an equilibrium system, the type of switching process described above is an 
activation process with a residence time described by the Arrhenius law. A simple model of 
such an activation process is the Langevin dynamics of a Hamiltonian system in a potential 
V. The noise level is then related to the temperature, and the Arrhenius law takes the 



form 22,35 37 



r;^^ ocexp(AV+_/a), (41) 

r:"' ocexp(AV_+/a). (42) 

Here, AV+_ and AV_+ are respectively the potential energy barrier as observed from the 
inhomogeneous and the homogeneous state. In a non-equilibrium context such as ours, there 



is no obvious equivalent of a potential, but the law given by Eqs. (41) and (42) is expected 



to hold on a fairly general basis, in the limit of small noise. This may be established from the 
instanton theory, or, from the Freidlin-Wentzell theory, which allows to compute V explicitly 
for a given model 38 , 39 . Our system does not fulfill the hypothesis of Freidlin-Wentzell 



theory, nevertheless, it is interesting to check if the law given by Eqs. (41) and (42) holds 



Our simulation data shown in Fig. 9 show that the dependence of r^*^* on a, as in Eqs. (41) 



and (42), holds also for our model, thereby suggesting that in the limit of low noise, the 



system behaves as one with transitions activated by a weak noise. 
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Figure 9. The plot shows as a function of 1/a the log of the mean residence time r 
in the two bistable states, namely, the inhomogeneous (m > 0) state and the homogeneous 
state (m w 0). The plot is based on data obtained from simulations with modes 1 — 3 
excited, whose amplitudes satisfy ci — C2 ~ = C(0)/6, where C(0) = 2T. Here, the 
kinetic temperature T = 0.53, while At = 10^^, iV — 5000. The straight line fits imply 
that rl^i ^ exp(l/Q;), in accordance with Eqs. (41) and (42). That the slopes of the two 



straight lines in the plot are different could be due to the fact that the height of the barrier 
is different when observed from the inhomogeneous and the homogeneous state. 



We conclude this section by describing briefly the algorithm to find P±{t^) to produce 
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Fig. [8| and t^'^ to produce Fig. 9j To this end, one has to identify from the time series data 
of the magnetization (see Fig. 7](a) for an example) the switching time instants between 
the two states. In the hmit of very small a, the distinction between two states should be 
obvious. However, we could not reach such a limit in our numerical simulations because the 
simulation time grows exponentially with 1/a (see Fig. [9]). For intermediate values of a, it is 
then a challenge to define precisely the two states. Indeed, as may be seen in Fig. [7](a), the 
data show strong fluctuations and hence, one needs to filter out "spurious" switching events 
and retain only the genuine ones. This may be done efficiently as we now discuss. 

We first obtain from the data a rough estimate of the mean of the magnetization when 
the system is in the two states, the homogeneous and the inhomogeneous one. Let us 
denote by m> and m< these estimates when the system is in the inhomogeneous and the 
homogeneous state, respectively. Let us define a "threshold" value of the magnetization rrith 
as the average of m> and m<; magnetization crossing this threshold to switch from one state 
to another however is not a precise enough criterion to define a switching event, as is obvious 
from Fig. [7j We thus resort to our algorithm which we now illustrate for the case when 
the system is in the homogeneous state; when the system is in the inhomogeneous state, 
the algorithm may be defined in a manner similar to the one below. In our algorithm, we 
identify a switching event as the one for which the following two conditions are satisfied, 
namely, (i) that the magnetization crosses the threshold to switch from the homogeneous to 
the inhomogeneous state, (ii) the magnetization after the switching reaches the value m> 
before reaching the value m<. When a switching event occurs, the switching time is defined 
as the time at which the magnetization crossed the threshold. This algorithm allows us to 
precisely define the switching times, from which we compute the switching time statistics 
P±{r±), and hence, the mean residence time r^*^*. 



6. Conclusions 



In this work, we considered long-range interacting systems driven by external stochastic 
fields, thereby leading to generic nonequilibrium stationary states. To study spatially 
homogeneous stationary states, we developed a kinetic theory approach by generalizing the 
known results for isolated long-range systems. Our theoretical approach is quite general, 
being applicable to any long-range inter-particle potential, space dimensions and boundary 
conditions. Our extensive numerical simulations on a paradigmatic model of long-range 
interacting systems demonstrated a very good agreement with the theory. Besides, our 
simulations for this representative case illustrated very interesting bistable behavior between 
homogeneous and inhomogeneous states, with a mean residence time that diverges as an 
exponential in the inverse of the strength of the external stochastic forces in the limit of low 
values of such forces. 

Let us note that another route to deriving the kinetic theory studied in this paper 
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is to adopt an approach similar to the one due to Khmontovich for isolated systems, 
by writing down the time evolution equation for the noise-averaged empirical measure 
p{p,q,t) = {^/N)J2iLi{^ilii^) ~ l)^{Pi{^) ~ P))- the resulting equation, the noise 
appears as a multiplicative term, which can be treated perturbatively, leading to the kinetic 
equation (32). 

This work leaves open some interesting issues, e.g., for technical simplicity, we assumed 
a homogeneous stationary state for the development of the kinetic theory. It would be 
of interest to generalize the theory to inhomogeneous states; in this regard, the method 



due to Heyvaerts reported recently may come of help 40 . Another issue is to study the 



dynamics of the kinetic equation (32), both analytically and numerically, which may unveil 
very interesting behaviors, such as limit cycles. One may also hope to develop a kinetic 
theory similar to the one analyzed here for related systems, for example, the point vortex 
model and the Euler equations in two-dimensional turbulence [Sl. 
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Appendix A. Condition of detailed balance for the dynamics ([s]) 

We prove here that the dynamics defined by the equations of motion ^ satisfies detailed 
balance if and only if = c for all fc, that is, if the stochastic forcing has a white spectrum 
in space. 

We start from the iV-particle Fokker-Planck equation ( 14 ) associated with the equations 
of motion It will be useful to rewrite it in the following way: 

^ - - E |-|A(x)/.(x)l + \ t ^IB,A^}M^)] , (A.1) 

1=1 *ii=i 

where Xi = qi for i = I, A^, Xi = Pi^N for i = {N + 1), ...,2N, and we use the notation 
X = {xi}. The drift vector Aj(x) is a function of the Xj's, and is given by 

A,{^)=Pi for t = l,...,N, (A.2) 

A.(x) = -ap^^j, - 1 V ~ ^^'^ for z = (AT + 1), 2N. (A.3) 

N ^ dqi_N 
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Similarly, the expression for the (symmetric) diffusion matrix Bi j is: 

B,j{x) = aCilqi.N - Qj-nD for i > N A j > N, (A.4) 
and Bij{x.) = otherwise. We moreover introduce the constants Ei = ±1, which denote the 
parity with respect to time inversion of the variables Xj, and the notation ex = {siXi}. 

Sect. 6.4) that the dynamics described 



It can be shown (see 22 , Sect. 5.3.5, or 25 



by a Fokker-Planck equation of the form (A.l) satisfies detailed balance if and only if the 

= l,...,2iV): 



following two conditions are satisfied 

B 



EiS jBi j (^^x) 



and 



2N 



e.A{e^)fU^) = -A,(x)/^(x) + ^ 



d 



1 dxj 



(A.5) 



(A.6) 



where /Xr(x) is the stationary solution of ( |A.1[ ). 



In our case, in which the drift and the diffusion terms are given by Eqs. (A. 2) and (A.4), 
respectively, the condition (A.5) is trivially satisfied. Our proof goes as follows: we solve 
formally Eq. (A.6), and show that f^{x.) is a stationary solution of Eq. (A.l) if and only if 
the non-vanishing part of Bij is proportional to the identity matrix. Then, it is simple to 
show that this implies that the spectrum of the forcing has to be white in space. 

Equation (A.6) for i = 1, ...,N is also trivially satisfied. On the other hand, for what 
(iV + 1),...,2A^, we have 

N 



concerns i 



dp. 



(A.7) 



where k = i — N. We introduce the N x N matrix C whose components are given by 
Cfcj(x) = C{\qk — qj\), and observe that, for generic values of the g^'s, C admits an inverse 
C~^. Integrating Eq. (A.7), we thus have 

N 

/^(x) = rf(gi,...,gjv)exp - ^ pk {C'^)^ .pj 

- k,j=l 

where d{qi, ...,qN) is an undetermined function. Inserting Eq. (A.8) into the Fokker-Planck 
equation ( |A.1[ ), imposing that it is a stationary solution, and with some calculations, we get 



(A.8) 



N 

E 



=1 >- 



NdH 



+ 



dH a/ 



N 



dqi dpi dqi dp^ 



0. 



Then, is a function of the Hamiltonian if, that is ffq{^) 



(A.9) 

iIj{H{x.)) for some function ip. 



On the other hand, because is given by the formula in Eq. (A.8), we can also deduce that 
ip is an exponential, and thus, that is Gaussian in the velocities. We conclude that 
(and hence, C) has to be independent of the g^'s and proportional to the identity. Finally, 
from the form of C(|gj — g^l) in Eq. (|4|, we see that this condition on C is satisfied if and 
only if the spectrum of the forcing is white in space. 
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Appendix B. Closure of the BBGKY hierarchy (17) 



We analyze here in detail the closure of the BBGKY hierarchy discussed in the text, in 
particular, the reasons for which the connected part of the two-particle correlation is of 
order a, while higher correlations are negligible at leading order in a, so that this closure is 
self-consistent. 

In the following, we expand the functions /2 and /a as 

f2iZi,Z2,t) = f{zi,t)f{z2,t) +g{Zi,Z2,t), (B.l) 

and 

fsizi, Z2, Z3, t) = f{zi,t)f{z2, t)f{z3, t) + f{zi, t)g{z2, ^3, t) 

+ fiz2, t)g{zi, zs, t) + f{z3, t)g{zi, Z2, t) + h{zi, Z2, z^, t),(B.2) 

and similarly, for other /^'s for s > 4. 

Now, let us write explicitly the first two equations of the BBGKY hierarchy (17). The 
first one, obtained from Eqs. (17) and ( B.l[ ), is 

where 



dt dq dp dq 



Qp2 Qp 



dgidpi v'{q-qi)g{z, Zi, t), (B.3) 



dqidpiv{q- qi)f{qi,pi,t) 



(B.4) 



dRSl) to get 



is the mean-field potential. For the second equation of the hierarchy, we use Eqs. ( B.2[ ) and 

dg ^ dg d<l>[f] ^ f{z2)dv{qi-q2)df ^ 1 dvjqi 
dqi dpi dqi N dqi dpi 

df 



dg{zi,Z2,t) 
dt 



(I2) dg 



dpi N dqi 



+ 



dpi 
df df 



dpi 
+ {1^2}, 



dviqi-q-i) ^ d a 

dzs^^giz2,Z3) + g^^lapig] + ^C{\q, - q2\)g^^g^^ 



d'g 



dpidp2 



dz- 



dvjqi - ga) dh 
dqi dpi - 



(B.5) 



where the symbol {1^2} stands for an expression obtained from the bracketed one on the 
right hand side by exchanging the subscripts 1 and 2. 

Let us analyze the order of magnitude of various terms in Eq. (B.5). First of all, we 
have / 



1, as it is normalized to unity. However, we do not know a priori the order of 

91-92) df 
dqi dpi 



magnitude of g and h. Thus, the order of magnitude of all but the terms /(^ ^^(^i g2) _g/ 



and fC(|gi-g2|) 



dpi dp2 



is unknown. In the continuum limit Na ^ 1, we have 



f{z2) dv{qi - ga) df 
N dqi dpi 



N 



a ,, df df 

<^ar^ -C{\qi - gal)- 



dpi dp2 



(B.6) 
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so that it is natural to guess that g 
have 

f{z2) dv{qi - ga) df 



~ a. Let us also observe that in the limit Na <^ 1, we 

df df 



N 



dqi 
1/N. 



> a 



a 



C(|gi-g2|)- 



(B.7) 



dpi N 2 dpi dp2 

In the limit Na ^ 1, the kinetic theory leads to the Lenard- 



so that we obtain g 
Balescu equation. 

Once we have established that g ^ a, one can write down the equation of the hierarchy 
for h and, with similar reasoning as above, one then finds that h is at least of order a/N <^ a 
(or, depending on whether a/N <^ or the reverse), so that the term J dz^ ^^^g^"'^^^ ^ 



is negligible in Eq. (B.5), as may be straightforwardly checked. The iterative procedure 



can be repeated at all orders of the hierarchy. Discarding three-particle and higher-order 



correlations is thus a self-consistent procedure. Moreover, note that in Eq. (B.5), some of 
the terms are of higher orders (a^, a/N,...) with respect to a, and thus, can be discarded. 
The final form of the second equation of the BBGKY hierarchy is thus 



dg{zi,z2) 
dt 



d~g d~g dm 

~P^ir~ + 7i — ^ — + 

oqi opi 



q2\) 



dqi 
dldl 
dpi dp2 



dpi 



dzWiqi - q^)g{z2,z-i) 

+ {1^2}. 



Note that g ^ a implies, see Eq. (B.3), that the mean-field effect of the stochastic forces 



gives a contribution at the same order to the two-particle correlation induced by them. 



Appendix C. Evolution of the kinetic energy for the dynamics ([s]) 

We derive here the evolution of the kinetic energy as obtained from the kinetic equation 



(32). Let us recall that the average kinetic energy density at time t in the continuous limit 



can be written as 



dpp^ f{p,t). 



id] 



The starting point to obtain its time evolution is to multiply the kinetic equation (32) by 



1^2 



IP", and then, to integrate over p. Neglecting for the moment the non-linear part of the 
diffusion coefficient, and integrating by parts, we get 



C(0) = 0, 



which gives 



{km - 



cm 



4 / ^ 4 ' 



(C.2) 



(C.3) 



The kinetic energy density in the stationary state is thus {n)^^ = C(0)/4. 



We now have to prove that the non-linear part of the diffusion coefficient (33) does 



not contribute to the time evolution of the kinetic energy. Such a result is expected and is 
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usually valid for coUisional terms (i.e., those terms in the kinetic equations which are given by 
two-particle correlation), for example, in the Boltzmann equation or in the Lenard-Balescu 



equation 28 . The contribution to k{t) from the non-linear part of the diffusion coefficient 



is a sum of terms proportional to 
1 

2 _ 

d 



T 



X 



dp 



dpp^ 



\e{k, A;p)p 



dpi h 

Pi-p 



dpi 



f'iPi,t) 



Pi-p |e(fc,/cpi)|2_ 



(C 



we will show that each of such terms vanishes independently. Indeed, integrating the last 
expression over p by parts, we get that 



T 



dpipf'{p,t) 



f'iPut) 



+ 



f'iPut) 



.(C.5) 



Pi — P l^ik, kp)\'^ Pi — p \e{k,kpi)\'^ 
Exchanging now the variables pi and p and the order of integration, we get that the above 
equation may be rewritten as 



T = dp dpipi f'{pi,t) 



f'{p,t) 



+ 



Pi — p \e{k^kpi)\^ Pi — P \£{k,kp)\'^ 
Summing up the last two equations, we therefore have 



T 



dp / dpif'ipi,t)fip,t) 



+ 



1 



_\e{k,kpi)\'^ \e{k,kp)\ 



. (C.6) 



(C.7) 



which vanishes on integrating by parts both with respect to pi and p. 



Appendix D. Proof that Eq. (32) admits non-Gaussian stationary distribution 
with Gaussian tails 

We prove here that for a general forcing spectra, the Gaussian distribution function in 



Eq. (37) is not a stationary solution of the kinetic equation (32), and that the tails of any 



stationary state are Gaussian. For the first point, we have to prove that the contribution to 



df/dt from the non-linear part of D[f]{p) in Eq. (33) is not vanishing. This result can be 
proven with an asymptotic expansion |12| for large momenta of the integrals which appear 
in the diffusion coefficient. Given any function g{p), we approximate integrals of the form 

'^P'^ (D.l) 



dpi- 

Pi-p 

by expanding in Taylor series. We get, for example. 



pi-p 



dp 



faiP 



Pi-P 



9 roo 



Pi 
P 



+ 



+ ... 



(D.2) 
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where, in the last equahty, we have taken into account the fact that the Gaussian distribution 

^ j with k odd do not contribute. In a similar way, we 



being even, the terms containing 
have 



\eik,kp)\' 



1 



and 



dpi 



1 



Pi — P l^ik, kpi 



2^3/2 



dpi 



pje- 



\e{k,kpi)\'^' 



(D.3) 



(D.4) 



where we have used the fact that \6{k,kp)\'^ is an even function of p. With these results, 
we can evaluate the non-linear part of the kinetic equation: for large pi, the non-linear 
contribution to df/dt is 



k=l 



1 + 



2/33/2 



vr 



dp 



p^e ^ 
e{k, A;p)p 



4/35/2 



TT 



(D.5) 



It can be shown that such a term is a non- vanishing function of pi. This completes the proof: 
for a generic forcing spectra, the stationary state, when exists, is not Gaussian. 

Using the same asymptotic expansion as before, it can be checked that the diffusion 
coefficient D[f]{p) converges to C(0)/2 for any distribution /. From this observation and 
Eq. (39), it follows that any stationary solution of the kinetic equation (32) has Gaussian 
tails. 
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